NeXT Education Software Sampler 1992 Fall/NeXT Education Software Sampler 1992 Fall.iso/Mathematics/Notebooks/SigProc2.0/Packages/SignalProcessing/Support/LatticeTheory.m

Text File | 1992-08-18 | 47.3 KB | 1,410 lines

(*  :Title:    Routines to implement concepts from Lattice Theory  *)

(*  :Authors:    Brian Evans, James McClellan  *)

(*
    :Summary:   Provide some of the functions which are common in lattice
        theory (such as computing distinct coset vectors) and the
        geometry of numbers (such as the Bezout identity and
        and generalized Euclid factors) that are not implemented
        in Mathematica.  These routines form the basis for
        multirate signal processing.
 *)

(*  :Context:    SignalProcessing`Support`LatticeTheory`  *)

(*  :PackageVersion:  2.7    *)

(*
    :Copyright:    Copyright 1989-1992 by Brian L. Evans
        Georgia Tech Research Corporation

    Permission to use, copy, modify, and distribute this software
    and its documentation for any purpose and without fee is
    hereby granted, provided that the above copyright notice
    appear in all copies and that both that copyright notice and
    this permission notice appear in supporting documentation,
    and that the name of the Georgia Tech Research Corporation,
    Georgia Tech, or Georgia Institute of Technology not be used
    in advertising or publicity pertaining to distribution of the
    software without specific, written prior permission.  Georgia
    Tech makes no representations about the suitability of this
    software for any purpose.  It is provided "as is" without
    express or implied warranty.
 *)

(*  :History:    multirate signal processing  *)

(*  :Keywords:    downsampling, sampling matrix, upsampling *)

(*
    :Source:         R. Bamberger, {The Directional Filter Bank:  A Multirate
        Filter Bank for the Directional Decompostion of Images},
        Georgia Tech Ph. D. Thesis, 1990.

             T. Kailath, {Linear Systems}, Prentice Hall, Englewood
        Cliffs, NJ, 1980.

             A. Kaufmann and A. Henry-Labordere, {Integer and
        Mixed Progamming: Theory and Applications}, Academic Press,
        New York, 1977.

             J. Kovacevic and M. Veterli, "Perfect Reconstruction
        Filter Banks with Rational Sampling Rates in One- and
        Two-dimensions," in {Proc. SPIE Conference on Visual
        Communications and Image Processing} (Philadelphia, PA),
        pp. 1258-1268, November 1989.

             C. MacDuffee, {The Theory of Matrices}, Springer-Verlag,
        Berlin, Germany, 1933.
 *)

(*  :Warning:    *)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*
    :Discussion:     The built-in integer division operation Quotient does
        not work properly under Mathematica 1.2.  Therefore, we have
        added a local routine called intdivide.

             These routines underlie the "algebra" for multidimensional
        multirate systems.  They are a combination of concepts from
        matrix theory and the geometry of numbers.

             The lone Mathematica lattice operation, LatticeReduce,
        is the Lenstra-Lenstra-Lovasz lattice reduction algorithm
        which makes the basis vectors of a lattice more orthogonal.
        Another way to perform this operation is to compute the
        Smith Form decomposition of the integer matrix and discard
        the right regular unimodular matrix factor
 *)

(*
    :Functions:    BezoutNumbers
        ColumnHermiteForm
        CommutableResamplingMatricesQ
        DiagonalMatrixQ
        DistinctCosetVectors
        GCLD
        GCRD
        IncList
        ImpulseTrain
        IntegerVectorQ
        LCLM
        LCRM
        NormalizeSamplingMatrix
        RegUnimodularResMatrixQ
        RelativelyPrime
        ReorderResampling
        ResamplingMatrix
        ResamplingMatrixMod
        RowHermiteForm
        SmithFormSameU
        SmithFormSameV
        SmithNormalForm
        SmithOrderedForm
        SmithReducedForm
        UpsampledSystem
 *)


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ]


(*  Only for the release of LatticeTheory as an independent package  *)

(*
If [ SameQ[ Head[Dialogue::usage], MessageName ],
     Dialogue::usage =
    "Dialogue is an option for certain routines to justify their answers. \
    Possible settings are False, True, or All, for no, partial, or \
    full justification, respectively." ]
 *)


(*  B E G I N     P A C K A G E  *)


BeginPackage[ "SignalProcessing`Support`LatticeTheory`",
          "SignalProcessing`Support`SupCode`" ]


(*  U S A G E     I N F O R M A T I O N  *)

BezoutNumbers::usage =
    "BezoutNumbers[a, b] gives the Bezout numbers of integers a and b. \
    That is, it returns {lambda, mu} so that lambda a + mu b == gcd(a,b). \
    In this case, BezoutNumbers simply calls ExtendedGCD. \
    For the collection of integers a1, a2, ..., an, \
    BezoutNumbers[a1, a2, ..., an] returns a set (list) of of numbers \
    k1, k2, ..., kn such that k1 a1 + k2 a2 + ... + kn an = d, \
    where d is the greatest common divisor of a1, a2, ..., an. \
    See also ExtendedGCD."

ColumnHermiteForm::usage =
    "ColumnHermiteForm[F] returns {X, H} where X is the m x m regular \
    unimodular matrix that puts the m x n integer matrix F into upper \
    triangular (column Hermite) form H (which is m x n). \
    That is, X . F = H."

CommutableResamplingMatricesQ::usage =
    "CommutableResamplingMatricesQ[downsamplingmat, upsamplingmat] \
    returns True if the two resampling operations are commutable."

DiagonalMatrixQ::usage =
    "DiagonalMatrixQ[matrix] returns True if the matrix is diagonal, \
    and gives False otherwise."

DistinctCosetVectors::usage =
    "DistinctCosetVectors[resmat] returns a sorted list of all of the \
    distinct coset vectors associated with the resampling matrix resmat. \
    DistinctCosetVectors[resmat, U^-1, Lambda, V^-1] finds all of the \
    distinct coset vectors of resmat from its Smith Form  U Lambda V. \
    This is carried by finding the distinct coset vectors of the diagonal \
    matrix Lambda, mapping them by U^-1, and then taking each vector \
    modulo the resampling matrix resmat. \
    The returned coset vectors are not sorted."

EuclidFactors::usage =
    "EuclidFactors[k, p, q] gives the factors a and b such that \
    p a + q b == k for relatively prime p and q. \
    When p and q are integers, this routine searches over a set of \
    min(|p|, |q|) integers. \
    When p and q are integer matrices, this routine searches over a \
    set of min(|det(p)|, |det(q)|) integer vectors and the vector a \
    or b corresponding to the integer matrix with the smallest \
    determinant will always be a distinct coset of that matrix. \
    See also RelativelyPrime."

GCLD::usage =
    "GCLD[a, b] finds a greatest common left divisor of resampling \
    matrices a and b which is unique to within the product of \
    a regular unimodular matrix on the right. \
    It forms F = {{a, b}, {0, 0}} and solves the equation F . X = H \
    where X is a square regular unimodular matrix and H is the row \
    Hermite (lower-triangular) form of the rectangular matrix F. \
    Of the four submatrices of H, only H11 is non-zero. \
    The GCLD of a and b is H11. \
    This routine supports a Dialogue option. \
    See also LCLM, RowHermiteForm, and ResamplingMatrix."

GCRD::usage =
    "GCRD[a, b] returns a greatest common right divisor of resampling \
    matrices a and b which is unique to within the product of \
    a regular unimodular matrix on the left. \
    It forms F = {{a, 0}, {b, 0}} and solves the equation X . F = H \
    where X is a square regular unimodular matrix and \
    H is the column Hermite (upper-triangular) form of \
    the rectangular matrix F. \
    Of the four submatrices of H, only H11 is non-zero. \
    H11 is the GCRD. \
    This routine supports a Dialogue option. \
    See also LCRM, ColumnHermiteForm, and ResamplingMatrix."

IncList::usage =
    "IncList[list, start, end] will increment the first element of list. \
    If this becomes greater than the first element of end, \
    then the next element of list will be incremented, and so forth. \
    For example, IncList[{9,0,0}, {0,0,0}, {10,10,10}] would return \
    {0,1,0} and IncList[{9,9,9}, {0,0,0}, {10,10,10}] would return \
    {10,10,10}. \
    This is useful when enumerating values over an \
    n-dimensional rectangular prism."

ImpulseTrain::usage =
    "ImpulseTrain[D, n] represents an infinite collection of \
    Kronecker impulse functions. \
    In one-dimension, D and n are integers. \
    The multidimensional separable form is when both D and n are \
    integer vectors of the same length. \
    In the non-separable case, D is a resampling matrix and \
    n is an integer vector. \
    This function gives 1 whenever D^-1 n is an integer (vector), and \
    gives 0 otherwise. \
    See also ResamplingMatrix and Impulse."

IntegerVectorQ::usage =
    "IntegerVectorQ[arg] returns True if arg is a vector whose \
    elements are integers."

LCLM::usage =
    "LCLM[a, b] returns a least common left multiple of resampling \
    matrices a and b which is unique to within the product of \
    a regular unimodular matrix on the left. \
    It forms F = {{a, 0}, {b, 0}} and solves the equation X . F = H \
    where X is a regular unimodular matrix and H is the column Hermite \
    (upper-triangular) form of the rectangular matrix F. \
    Of the four submatrices of H, only H11 is non-zero. \
    The LCLM of a and b is (X21 . a) = -(X22 . b). \
    This routine supports a Dialogue option. \
    See also GCRD, ColumnHermiteForm, and ResamplingMatrix."

LCRM::usage =
    "LCRM[a, b] finds a least common right multiple of resampling \
    matrices a and b which is unique to within the product of \
    a regular unimodular matrix on the right. \
    It forms F = {{a, b}, {0, 0}} and solves the equation F . X = H \
    where X is a regular unimodular matrix and H is the row Hermite \
    (lower-triangular) form of the rectangular matrix F. \
    Of the four submatrices of H, only H11 is non-zero. \
    The LCRM of a and b is (a . X12) = -(b . X22). \
    This routine supports a Dialogue option. \
    See also GCLD, RowHermiteForm, and ResamplingMatrix."

NormalizeSamplingMatrix::usage =
    "NormalizeSamplingMatrix[sampmat] will decompose the sampling matrix \
    sampmat into a diagonal matrix D and a normalized sampling matrix V. \
    The results are returned as { D, V }."

ReorderResampling::usage =
    "ReorderResampling[resmat, vars, allvars] rewrites \
    the resampling matrix resmat and the variables vars associated \
    with it so that the resampling is carried out by the order \
    of variables specified by allvars. \
    Therefore, vars must be a subset of allvars."

RegUnimodularResMatrixQ::usage =
    "RegUnimodularResMatrixQ[mat] gives Ture if mat is a resampling \
    matrix with a determinant of +1 or -1, and gives False otherwise. \
    See also ResamplingMatrix."

RelativelyPrime::usage =
    "RelativelyPrime[a, b] returns True if a and b are relatively prime, \
    and gives False otherwise. \
    The data types of a and b must be the same. \
    Possible data types are integers, polynomials, and \
    square integer matrices (a.k.a. resampling matrices). \
    For non-diagonal resampling matrices, a third argument of \
    Right or Left is necessary to direct the matrix factorization \
    to occur on the right (using GCLD) or left (using GCRD). \
    See also GCD, PolynomialGCD, and ResamplingMatrix."

ResamplingMatrix::usage =
    "ResamplingMatrix[arg] returns True if arg is a square matrix \
    whose elements are integers and whose determinant is non-zero. \
    It returns False if arg is not a square matrix. \
    It also returns False if arg is a square matrix with one element \
    being a number that is not an integer. \
    Otherwise, it returns the conditions for which the square matrix \
    would be a resampling matrix."

ResamplingMatrixMod::usage =
    "ResamplingMatrixMod[n, N] implements the integer vector n \
    modulo a sampling matrix N which is defined as \
    n - N . Floor[ N^-1 . n ]. \
    The operation will be carried out if and only if n is an integer \
    vector, N is (or could be) a sampling matrix, and \
    N . n is a valid operation."

RowHermiteForm::usage =
    "RowHermiteForm[F] returns {X, H} where X is the n x n regular \
    unimodular matrix that puts the m x n integer matrix F into \
    lower triangular (row Hermite) form H (which is m x n). \
    That is, F . X = H."

SmithFormSameU::usage =
    "SmithFormSameU[smithfun, m1, m2] computes the Smith Form of the \
    integer matrix m1 so that m1 = u1 d1 v1. \
    The routine then uses u2 = u1 to decompose m2. \
    The first argument specifies which decomposition routine should \
    be used to decompose m1 (e.g., SmithNormalForm). \
    The routine returns {cond, {u1^-1, d1, v1^-1}, {u2^-1, d2, v2^-1}}, \
    where cond is True of False depending on whether or not matrices \
    m1 and m2 can be decomposed using the same U."

SmithFormSameV::usage =
    "SmithFormSameV[smithfun, m1, m2] computes the Smith Form of the \
    integer matrix m1 so that m1 = u1 d1 v1. \
    The routine then uses v2 = v1 to decompose m2. \
    The first argument specifies which decomposition routine should \
    be used to decompose m1 (e.g., SmithNormalForm). \
    The routine returns {cond, {u1^-1, d1, v1^-1}, {u2^-1, d2, v2^-1}}, \
    where cond is True of False depending on whether or not matrices \
    m1 and m2 can be decomposed using the same V."

SmithNormalForm::usage =
    "SmithNormalForm[A] returns a list of three matrices \
    U^-1, D, and V^-1 such that the m x n resampling matrix A \
    equals U D V where U (m x m) and V (n x n) are integer matrices \
    with determinants of +1 or -1 and D is a diagonal matrix \
    of integer components such that |det(D)| = |det(A)|. \
    The algorithm follows Kaufmann's {Integer and Mixed Programming}. \
    Note that this factorization is not unique. \
    SmithNormalForm supports a Dialogue option. \
    If Dialogue is True or All, SmithNormalForm will show immediate steps. \
    See also SmithOrderedForm and SmithReducedForm."

SmithOrderedForm::usage =
    "SmithOrderedForm[A] returns a list of three matrices \
    U, D, and V such that the m x n resampling matrix A \
    equals U D V where U (m x m) and V (n x n) are integer matrices \
    with a determinant of 1 and D is a diagonal matrix \
    of integer components such that det(D) = det(A). \
    The components along the diagonal of D are sorted by absolute value. \
        The diagonal elements are positive except for possibly the \
    last one whose sign will be the sign of det(A). \
    Note that this factorization is unique if none of the diagonal \
    elements are equal in absolute value. \
    Also note that the starting point is the result of the Smith \
    Normal Form decomposition. \
    Like SmithNormalForm, SmithOrderedForm supports a Dialogue option. \
    If True or All, SmithNormalForm will show immediate steps. \
    See also SmithNormalForm and SmithReducedForm."

SmithReducedForm::usage =
    "SmithReducedForm[A] returns a list of three matrices \
    U, D, and V such that the m x n resampling matrix A \
    equals U D V where U (m x m) and V (n x n) are integer matrices \
    with a determinant of +1 or -1 and D is a diagonal matrix \
    of positive integers such that |det(D)| = |det(A)|. \
    Each component down the diagonal is a factor of the next one. \
    As a result, the diagonal elements are sorted down the diagonal. \
    The factorization is not unique. \
    Note that the starting point for this algorithm is the Smith \
    Ordered Form. \
    Like SmithOrderedForm, SmithReducedForm supports a Dialogue option. \
    If True or All, SmithNormalForm will show immediate steps. \
    See also SmithNormalForm and SmithOrderedForm."

UpsampledSystem::usage =
    "UpsampledSystem[expr, vars] returns a resampling matrix and \
    a set of linearly coupled variables taken from the master list vars. \
    For example, UpsampledSystem[Cos[w1 + w2] Sin[w1 - w2], {w1, w2} ] \
    will return { {{1,1},{1,-1}}, {w1,w2} }. \
    In this case, the dot product of the resampling matrix and the \
    variable list give {w1 + w2, w1 - w2}. \
    A third optional argument would be a function to extract coefficients \
    that takes an expression, a list of variables, and a value to return \
    on failure."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin[ "`Private`" ]


(*  M E S S A G E S  *)

BezoutNumbers::error =
    "BezoutNumbers[``, ``] with gcd = `` could not be found."

CommutableResamplingMatricesQ::invalid =
    "At least one of the arguments is not a resampling matrix."

DistinctCosetVectors::notresmat =
    "The lone argument is not a resampling matrix: ``"
DistinctCosetVectors::zerodet = "Resampling matrix has zero determinant."

EuclidFactors::notrelprime = "`` and `` are not relatively prime ``."

LCLM::badgcd = "The (``,``) position of F is not the gcd of column i"

SmithNormalForm::invalid =
    "The determinants of U^-1 and V^-1 should both be one."


(*  S U P P O R T I N G     R O U T I N E S  *)


(* identitymatrix *)
identitymatrix[{k_Integer, k_Integer}] := IdentityMatrix[k]
identitymatrix[{m_Integer, n_Integer}] :=
    Block [    {i, imat, zerovector},
        zerovector = Table[0, {n}];
        imat = Table[ zerovector, {m} ];
        For [ i = 1, i <= m, i++,
              imat[[i,i]] = 1 ];
        imat ] /;
    ( m != n )

(* intdivide --  kludge for Quotient primitive which has a problem  *)
(*         under Mma 1.2 when the second argument is negative *)
If [ TrueQ[ $VersionNumber >= 2.0 ],
     intdivide[n_, m_] := Quotient[n, m],
     intdivide[n_, m_] := Sign[m] Quotient[n, m Sign[m]] ]

(* mergerows *)
mergerows[ {x__} ] := Join[x]

(* matrix4to1 --  combine four component matrices into 1 *)
matrix4to1[a11_, a12_, a21_, a22_] :=
    Map[ mergerows, Transpose[{Join[a11, a21], Join[a12, a22]}] ]

(* matrix1to4 --  split one matrix into its four m x n component matrices *)
matrix1to4[mat_] := matrix1to4[ mat, Dimensions[mat] ]
matrix1to4[mat_, {m_, n_}] :=
    Block [    {a11, a12, a21, a22, i, j},
        a11 = Table[ mat[[i,j]], {i, 1, m}, {j, 1, n} ];
        a12 = Table[ mat[[i,j]], {i, 1, m}, {j, n+1, 2 n} ];
        a21 = Table[ mat[[i,j]], {i, m+1, 2 m}, {j, 1, n} ];
        a22 = Table[ mat[[i,j]], {i, m+1, 2 m}, {j, n+1, 2 n} ];
        {{a11, a12}, {a21, a22}} ]

(* myabs *)
myabs[x_List] := Min[Abs[Select[x, nonzero]]]

(*  mymod  *)
SetAttributes[mymod, {Listable}]
mymod[x_] := Mod[x, 1]

(* nonzero *)
nonzero[x_] := ( x != 0 )

(*  onefun  *)
onefun[x__] := 1

(*  toList  *)
toList[] := {}
toList[arg_List] := arg
toList[arg_] := List[arg] /; ! SameQ[Head[arg], List]
toList[arg1_, args__] := List[arg1, args]

(*  vectorgcd   *)
vectorgcd[ x_ ] := Apply[GCD, x]


(* FindMinMatrixElement *)
FindMinMatrixElement[mat_List] :=
    Block [ {item},
        item = myabs[ Flatten[mat] ];
        Join[ Position[ mat, item ], Position[ mat, -item ] ] ]

(* LCLMandGCRDintroduction *)
LCLMandGCRDintroduction[] :=
    Block [ {matrixborder},
        matrixborder = MatrixForm[ {"|", "|"} ];
        Print[ "  This routine builds the matrix X such that" ];
        Print[ "     ", matrixborder,
               MatrixForm[{{"X11", "X12"}, {"X21", "X22"}}],
               matrixborder, "  ", matrixborder,
               MatrixForm[{{"A", 0}, {"B", 0}}],
               matrixborder, Superscript[" = "], matrixborder,
               MatrixForm[{{"H", 0}, {0, 0}}], matrixborder ];
        Print[ "where A and B are the two resampling matrices." ];
        Print[ "X is composed of a product of unimodular" ];
        Print[ "matrices.  The least common left multiple is:" ];
        Print[ "M = X21 A = -X22 B and the greatest common" ];
        Print[ "right divisor of A and B is H." ] ]

(* LCRMandGCLDintroduction *)
LCRMandGCLDintroduction[] :=
    Block [ {matrixborder},
        matrixborder = MatrixForm[ {"|", "|"} ];
        Print[ "  This routine builds the matrix X such that" ];
        Print[ "     ", matrixborder,
               MatrixForm[{{"A", "B"}, {0, 0}}],
               matrixborder, "  ", matrixborder,
               MatrixForm[{{"X11", "X12"}, {"X21", "X22"}}],
               matrixborder, Superscript[" = "], matrixborder,
               MatrixForm[{{"H", 0}, {0, 0}}], matrixborder ];
        Print[ "where A and B are the two resampling matrices." ];
        Print[ "X is composed of a product of unimodular" ];
        Print[ "matrices.  The least common right multiple is:" ];
        Print[ "M = A X12 = -B X22 and the greatest common" ];
        Print[ "left divisor of A and B is H.  By transposing" ];
        Print[ "both sides of the equation, the solution steps" ];
        Print[ "are the same as carried out by the LCLM and" ];
        Print[ "GCRD routines.  This how the LCRM and GCLD" ];
        Print[ "routines work." ] ]

(* TransitionMatrix *)
TransitionMatrix[r1_, r2_, dim_] :=
        TransitionMatrix[r1, r2, dim, IdentityMatrix[dim]]
TransitionMatrix[r1_, r2_, dim_, mat_] :=
    Block [    {resmat, temp},
        resmat = mat;
        temp = resmat[[r1]];
        resmat[[r1]] = resmat[[r2]];
        resmat[[r2]] = temp;
        resmat ]


(* SmithFormA --  returns non-zero minimum value, left transition *)
(*   matrix, and right transition matrix of the current matrix.   *)
SmithFormA[ curmat_, dim_, m_, n_, identm_, identn_ ] :=
    Block [    {col, i, minvalue, row, workmat},
        workmat = Table[ Take[ curmat[[i]], {1 + dim, n} ],
                 {i, 1 + dim, m} ];
        {row, col} = First[FindMinMatrixElement[workmat]];
        { curmat[[row + dim]][[col + dim]],
          TransitionMatrix[dim + 1, row + dim, m, identm],
          TransitionMatrix[dim + 1, col + dim, n, identn] } ]

(* SmithFormB --  returns the left and right subtraction matrices *)
SmithFormB[ curmat_, dim_, minvalue_, m_, n_, identm_, identn_ ] :=
    Block [    {a, i, j, submatleft, submatright},

        submatleft = identm;
        For [ i = dim + 2, i <= m, i++,
              a = curmat[[i, 1 + dim]];
              If [ ! SameQ[a, 0],
                   submatleft[[i, 1+dim]] = intdivide[-a, minvalue]]];

        submatright = identn;
        For [ j = dim + 2, j <= n, j++,
              a = curmat[[1 + dim, j]];
              If [ ! SameQ[a, 0],
                submatright[[1+dim, j]] = intdivide[-a, minvalue]]];

        { submatleft, submatright } ]

(* SmithFormEndTest --  end if the elements in the first column  *)
(*   and first row of the working matrix except (1,1) are zero.  *)
SmithFormEndTest[curmat_, dim_, m_, n_] :=
    Block [ {firstrowcollist},

        firstrowcollist = Join[ Take[ curmat [[dim+1]], {dim+2, m} ],
                    Take[ Transpose[curmat] [[dim + 1]],
                          {dim + 2, n} ] ];

        SameQ[First[firstrowcollist], 0] &&
          Apply[Equal, firstrowcollist] ]

(* SmithFormTrace *)
SmithFormTrace[ dialogue_, str_, mat_ ] :=
    If [ dialogue || SameQ[dialogue, All],
         Print[ ]; Print[ str ]; Print[ MatrixForm[mat] ] ]

SmithFormTrace[ dialogue_, str_ ] :=
    If [ dialogue || SameQ[dialogue, All], Print[ ]; Print[ str ] ]

SmithFormTrace2[ dialogue_, str_, mat_ ] :=
    If [ SameQ[dialogue, All],
         Print[ ];
         Print[ "  ", str, " subtraction matrix:" ];
         Print[ "  ", MatrixForm[mat] ] ]


(*  E X P O R T E D     R O U T I N E S  *)

(*  BezoutNumbers --  only examine one period of possible values  *)
(*    to find more than two Bezout numbers, the key to do it by   *)
(*    induction:  generate the Bezout numbers for the first two;  *)
(*    then do it for the GCD of the first two no. and the third;  *)
(*    adjust the first two coefficients accordingly...            *)
BezoutNumbers[a_Integer, b_Integer] := ExtendedGCD[a, b] [[2]]

BezoutNumbers[a_Integer, b_Integer, rest__] :=
    Block [ {bezlist, intlist, lastgcd, j, n, newgcd, term},
        intlist = {a, b, rest};
        n = Length[intlist];

        {lastgcd, bezlist} = ExtendedGCD[a, b];
        For [ j = 3, j <= n, j++,
              term = intlist[[j]];
              newgcd = GCD[lastgcd, term];
              {p, q} = bezoutwithgcd[lastgcd, term, newgcd];
              bezlist *= p;
              AppendTo[bezlist, q];
              lastgcd = newgcd ];
        bezlist ] /;
    IntegerVectorQ[ {a, b, rest} ]

bezoutwithgcd[0, b_Integer, gcd_Integer] := {0, Sign[b]}
bezoutwithgcd[a_Integer, 0, gcd_Integer] := {Sign[a], 0}

bezoutwithgcd[a_Integer, b_Integer, gcd_Integer] :=
    Block [ {anorm, bnorm, lambda, mu, muvalue, mumax},

        anorm = a / gcd;
        bnorm = b / gcd;
        mumax = Abs[anorm];
        For [ mu = 0, mu < mumax, mu++,
              lambda = ( 1 - mu bnorm ) / anorm;
              If [ IntegerQ[lambda], muvalue = mu; Break[] ] ];

        If [ ! IntegerQ[lambda],
             Message[ BezoutNumbers::error, a, b, gcd ] ];

        { lambda, muvalue } ] /;
    ( a <= b )

bezoutwithgcd[a_Integer, b_Integer, gcd_Integer] :=
    Block [ {anorm, bnorm, lambda, lambdavalue, mu},

        anorm = a / gcd;
        bnorm = b / gcd;
        lambdamax = Abs[bnorm];
        For [ lambda = 0, lambda < lambdamax, lambda++,
              mu = ( 1 - lambda anorm ) / bnorm;
              If [ IntegerQ[mu], lambdavalue = lambda; Break[] ] ];

        If [ ! IntegerQ[mu],
             Message[ BezoutNumbers::error, a, b, gcd ] ];

        { lambdavalue, mu } ] /;
    ( a > b )

(* ColumnHermiteForm *)
Options[ColumnHermiteForm] := { Dialogue -> False }

HermiteEndTest[{}] := True
HermiteEndTest[ithcol_] := SameQ[First[ithcol], 0] && Apply[Equal, ithcol]

ColumnHermiteForm[ mat_, op___ ] :=
    Block [    {allflag, col, dialflag, dialogue, h, hji, htrans, i,
         identm, ithcol, j, m, n, oplist, pivot, r, temp, u, x},

        oplist = toList[op] ~Join~ Options[ColumnHermiteForm];
        dialogue = Replace[Dialogue, oplist];
        allflag = SameQ[dialogue, All];
        dialflag = dialogue || allflag;

        {m, n} = Dimensions[mat];
        r = Min[m, n];
        x = identm = IdentityMatrix[m];        (* unimodular *)
        h = mat;                (* Hermite matrix *)

        If [ dialflag,
             Print[ "Putting the ", m, " x ", n, " matrix below into" ];
             Print[ "upper triangular (Hermite) form:" ];
             Print[ MatrixForm[h] ] ];

        For [ i = 1, i <= r, i++,

              If [ allflag, Print[ "Iteration #", i ] ];

              htrans = Transpose[h];
              ithcol = Take[ htrans[[i]], {i+1, m} ];
              While [ ! HermiteEndTest[ithcol],

                  (*   Move smallest element in column i *)
                  (* that are below (i,i) to (i,i).      *)
                  PrependTo[ ithcol, h[[i,i]] ];
                  {col} = First[ FindMinMatrixElement[ ithcol ] ];
                  col += i - 1;
                  If [ i != col,
                       temp = h[[i]];
                       h[[i]] = h[[col]];
                       h[[col]] = temp;
                       temp = x[[i]];
                       x[[i]] = x[[col]];
                       x[[col]] = temp ];
                  If [ allflag,
                   Print[ "Hermite Form:" ];
                   Print[ MatrixForm[h] ] ];

                  (* reduce elements in column i below (i,i) *)
                  pivot = h[[i,i]];
                  u = identm;
                  Clear[j];
                  For [ j = i+1, j <= m, j++,
                    u[[j,i]] = intdivide[-h[[j,i]], pivot] ];
                  h = u . h;
                  x = u . x;

                  If [ allflag,
                   Print[ "Hermite Form:" ];
                   Print[ MatrixForm[h] ] ];

                      ithcol = Take[Transpose[h][[i]], {i+1, m}] ] ];

        {x, h} ] /;
    MatrixQ[mat] && IntegerVectorQ[ Flatten[mat] ]

(*  CommutableResamplingMatricesQ, adapted from [Kovacevic & Vetterli, p. 4] *)
(*  A cascade of a downsampling by D1 and upsampling by D2 is commutable if  *)
(*  1. The matrix product is commutable, i.e.   D1 . D2 == D2 . D1         *)
(*  2. The sets A and B are equivalent where                     *)
(*                                         *)
(*                                      t -1                     *)
(*                  A = { exp(-2 pi j (D )   k), k in LatD1 }             *)
(*                                      1                     *)
(*                                         *)
(*                                     t   t -1                     *)
(*                  B = { exp(-2 pi j D  (D )   n), n in LatD1 }         *)
(*                                     2   1                     *)
(*                                         *)
(*  The points in the lattice associated with D1 (LatD1) are none other      *)
(*  than the set of fundamental cosets.                         *)

Options[ CommutableResamplingMatricesQ ] := { Dialogue -> False }

CommutableResamplingMatricesQ[ D1_, D2_, op___ ] :=
    Block [    {anglesofA, anglesofB, cosets, dialogue, D1tinv, vectors},

        (*  check arguments  *)
        If [ ! ( ResamplingMatrix[D1] && ResamplingMatrix[D2] ),
             Return[Message[CommutableResamplingMatricesQ::invalid]] ];

        (*  check if dialogue is enabled  *)
        dialogue = Replace[ Dialogue,
                    toList[op] ~Join~
                    Options[CommutableResamplingMatricesQ] ];

        (*  Condition #1  *)
        If [ ! SameQ[D1 . D2, D2 . D1], Return[False] ];

        (*  Condition #2  *)
        cosets = DistinctCosetVectors[D1];
        D1tinv = Inverse[ Transpose[D1] ];
        vectors = D1tinv . Transpose[cosets];
        anglesofA = Union[Transpose[mymod[vectors]]];
        anglesofB = Union[Transpose[mymod[Transpose[D2] . vectors]]];

        If [ dialogue,
             Print["angles associated with the downsampling matrix:"];
             Print[2 Pi anglesofA];
             Print["angles associated with both resampling matrices:"];
             Print[2 Pi anglesofB] ];

        SameQ[anglesofA, anglesofB] ]

(*  DiagonalMatrixQ  *)
DiagonalMatrixQ[mat_] := False  /; ! MatrixQ[mat]
DiagonalMatrixQ[mat_] :=
        Block [ {cond = False, diag, dims, i},
        dims = Dimensions[mat];
        If [ dims[[1]] == dims[[2]],
                     diag = Table[ mat[[i,i]], { i, 1, dims[[1]] } ];
             cond = SameQ[mat, DiagonalMatrix[diag]] ];
        cond ]

(*  DistinctCosetVectors--  resmat is a dim x dim matrix  *)
DistinctCosetVectors[ resmat_?ResamplingMatrix ] :=
    Apply[DistinctCosetVectors, Prepend[SmithNormalForm[resmat], resmat]]

DistinctCosetVectors[ resmat_, uinv_, lambda_, vinv_ ] :=
    Block [    {absdet, cosets, curvertex, diagvector, dim,
         invresmat, lowervertex, u, uppervertex, zerovector},

        absdet = Abs[Det[resmat]];
        If [ SameQ[absdet, 0],
             Message[ DistinctCosetVectors::zerodet ];
             Return[] ];

        dim = Apply[Min, Dimensions[resmat]];
        zerovector = Apply[Table, {0, {dim}}];
        cosets = { zerovector };

        If [ SameQ[absdet, 1], Return[cosets] ];

        invresmat = Inverse[resmat];
        u = Inverse[uinv];

        diagvector = Table[ lambda[[i,i]], {i, 1, dim} ];
        uppervertex = Abs[ diagvector ];
        lowervertex = curvertex = zerovector;

        (*  enumerate all points in the rectangular prism      *)
        (*    delimited by the column vectors of the diagonal     *)
        (*    matrix lambda, map them ny U^-1, and then modulo    *)
        (*    each one with the original resampling matrix resmat *)
        cosets = Table[ curvertex = IncList[ curvertex,
                             lowervertex,
                             uppervertex ];
                ResamplingMatrixMod[ u . curvertex,
                               resmat, invresmat ],
                {absdet - 1} ];

        Prepend[cosets, zerovector] ]

(*  EuclidFactors--  p a + q b == k  with p, q, and k given w/ GCD(p,q) == 1 *)
(*  1.  integer case  *)
EuclidFactors[0, p_Integer, q_Integer] := {0, 0}
EuclidFactors[k_Integer, p_Integer, q_Integer] :=
    If [ p < q,
         Reverse[ integerEuclidFactors[k, q, p] ],
         integerEuclidFactors[k, p, q] ]

integerEuclidFactors[k_, p_, q_] :=
    Block [    {a, b, continue, kmod, nump, numq, periods, pxq},
        pxq = p q;
        periods = intdivide[k, pxq];
        kmod = Mod[k, pxq];
        nump = intdivide[periods, 2];
        numq = nump + Mod[periods, 2];
        continue = True;
        a = 0;
        While [ continue,
            b = ( kmod - p a ) / q;
            If [ IntegerQ[b], continue = False, a++ ] ];
        {a + q nump, b + p numq} ] /;
    GCD[p, q] == 1

integerEuclidFactors[k_, p_, q_] :=
    Message[EuclidFactors::notrelprime, p, q, "integers"]

(*  2. Integer matrix case  *)
EuclidFactors[k_?IntegerVectorQ, p_?ResamplingMatrix, q_?ResamplingMatrix] :=
    Block [    {detp, detq},
        detp = Det[p];
        detq = Det[q];
        If [ Abs[detp] > Abs[detq],
             Reverse[ vectorEuclidFactors[k, q, p, detq, detp] ],
             vectorEuclidFactors[k, p, q, detp, detq] ] ]

vectorEuclidFactors[k_, p_, q_, detp_, detq_] :=
    Block [    {a, b, cosets, i, numcosets, pinv, pinvDOTk, pinvDOTq},

        (* We can rewrite  k = p . a  +  q . b   as               *)
        (*             -1             -1                          *)
        (*            p   . k  -  p  . q . b  =  a             *)
        pinv = Inverse[p];
        pinvDOTk = pinv . k;
        pinvDOTq = pinv . q;

        (* Now, we want to enumerate all possible values of b     *)
        (* and find the a that makes the above equation give      *)
        (* an integer vector.  We have mapped k into the          *)
        (* fundamental paralellopiped so all we have to do is     *)
        (* enumerate the cosets of q.                             *)
        cosets = DistinctCosetVectors[q];
        numcosets = Length[cosets];
          For [ i = 1, i <= numcosets, i++,
                b = cosets[[i]];
                a = pinvDOTk - pinvDOTq . b;
                If [ IntegerVectorQ[a], Break[] ] ];
        {a, b} ] /;
    RelativelyPrime[p, q, Right]

vectorEuclidFactors[k_, p_, q_, detp_, detq_] :=
    Message[ EuclidFactors::notrelprime, p, q,
         "(on the right) integer matrices" ]


(*  GCLD -- a and b are m x n; X is 2n x 2n; H and F are 2m x 2n     *)
Options[GCLD] := { Dialogue -> False }
GCLD[ a_, b_, op___ ] :=
    Block [ {allflag, dialflag, dialogue, F, H, H11, H12, H21, H22,
         m, n, oplist, X, zeromxn, zerovector},

        oplist = toList[op] ~Join~ Options[GCRD];
        dialogue = Replace[Dialogue, oplist];
        allflag = SameQ[dialogue, All];
        dialflag = dialogue || allflag;
        If [ dialflag, LCRMandGCLDintroduction[] ];

        {m, n} = Dimensions[a];
        zerovector = Table[0, {n}];
        zeromxn = Table[zerovector, {m}];
        F = matrix4to1[a, b, zeromxn, zeromxn];

        {X, H} = RowHermiteForm[F, oplist];
        {{H11, H12}, {H21, H22}} = matrix1to4[H, {m, n}];

        If [ allflag,
             Print[ "Hermite Form (right side):" ];
             Print[ MatrixForm[H] ];
             If [ ! SameQ[F . X, H],
              Print[ "Error:  F . X != H" ] ] ];

        H11 ] /;
    Dimensions[a] == Dimensions[b] && ResamplingMatrix[a] &&
      ResamplingMatrix[b]

(* GCRD -- a and b are m x n, X is 2m x 2m, and H is 2m x 2n *)
Options[GCRD] := { Dialogue -> False }
GCRD[ a_, b_, op___ ] :=
    Block [ {allflag, dialflag, dialogue, F, H, H11, H12, H21, H22,
         m, n, oplist, X, zeromxn, zerovector},

        oplist = toList[op] ~Join~ Options[GCRD];
        dialogue = Replace[Dialogue, oplist];
        allflag = SameQ[dialogue, All];
        dialflag = dialogue || allflag;
        If [ dialflag, LCLMandGCRDintroduction[] ];

        {m, n} = Dimensions[a];
        zerovector = Table[0, {n}];
        zeromxn = Table[zerovector, {m}];
        F = matrix4to1[a, zeromxn, b, zeromxn];

        {X, H} = ColumnHermiteForm[F, oplist];
        {{H11, H12}, {H21, H22}} = matrix1to4[H, {m, n}];

        If [ allflag,
             Print[ "Hermite Form (right side):" ];
             Print[ MatrixForm[H] ];
             If [ ! SameQ[X . F, H],
              Print[ "Error:  X . F != H" ] ] ];

        H11 ] /;
    Dimensions[a] == Dimensions[b] && ResamplingMatrix[a] &&
      ResamplingMatrix[b]

(*  IncList  *)
IncList[list_, start_, end_] :=
    Block [    {incflag, newlist},
        newdigit[d_, s_, e_] :=
            Which [ incflag && SameQ[d + 1, e],
                  s,
                incflag,
                  incflag = False;
                  d + 1,
                True,
                  d ];
        incflag = True;
        newlist = newdigit[list, start, end];
        If [ SameQ[newlist, start], end, newlist ] ]

SetAttributes[newdigit, Listable]

(*  ImpulseTrain  *)
ImpulseTrain[D_Integer, n_Integer] := 1        /; IntegerQ[n / D]
ImpulseTrain[D_Integer, n_Integer] := 0

ImpulseTrain[D_?ResamplingMatrix, n_?IntegerVectorQ] := 1    /;
    IntegerVectorQ[ Inverse[D] . n ]
ImpulseTrain[D_?ResamplingMatrix, n_?IntegerVectorQ] := 0

ImpulseTrain[{D_}, {n_}] := ImpulseTrain[D, n]
ImpulseTrain[D_List, n_List] :=
    Apply[Times, Map[ itrain, Transpose[ { D, n } ] ] ] /;
    VectorQ[D] && VectorQ[n] && Length[D] == Length[n]

itrain[{D_, n_}] := ImpulseTrain[D, n]

(*  IntegerVectorQ  *)
IntegerVectorQ[arg_] := VectorQ[arg] && Apply[And, Map[IntegerQ, arg]]

(*  LCLM -- a and b are m x n, X is 2m x 2m, and H is 2m x 2n  *)
Options[LCLM] := { Dialogue -> False }
LCLM[ a_, b_, op___ ] :=
    Block [ {allflag, dialflag, dialogue, F, H, m, n, oplist, X,
         X11, X12, X21, X22, zeromxn, zerovector},

        oplist = toList[op] ~Join~ Options[LCLM];
        dialogue = Replace[Dialogue, oplist];
        allflag = SameQ[dialogue, All];
        dialflag = dialogue || allflag;
        If [ dialflag, LCLMandGCRDintroduction[] ];

        {m, n} = Dimensions[a];
        zerovector = Table[0, {n}];
        zeromxn = Table[zerovector, {m}];
        F = matrix4to1[a, zeromxn, b, zeromxn];

        {X, H} = ColumnHermiteForm[F, oplist];
        {{X11, X12}, {X21, X22}} = matrix1to4[X, {m, m}];

        If [ allflag,
             Print[ "Hermite Form (right side):" ];
             Print[ MatrixForm[H] ];
             If [ ! SameQ[X . F, H],
              Print[ "Error:  X . F != H" ] ] ];

        X21 . a ] /;
    Dimensions[a] == Dimensions[b] && ResamplingMatrix[a] &&
      ResamplingMatrix[b]

(*  LCRM -- a and b are m x n; X is 2n x 2n; H and F are 2m x 2n     *)
Options[LCRM] := { Dialogue -> False }
LCRM[ a_, b_, op___ ] :=
    Block [ {allflag, dialflag, dialogue, F, H, m, n, oplist, X,
         X11, X12, X21, X22, zeromxn, zerovector},

        oplist = toList[op] ~Join~ Options[GCRD];
        dialogue = Replace[Dialogue, oplist];
        allflag = SameQ[dialogue, All];
        dialflag = dialogue || allflag;
        If [ dialflag, LCRMandGCLDintroduction[] ];

        {m, n} = Dimensions[a];
        zerovector = Table[0, {n}];
        zeromxn = Table[zerovector, {m}];
        F = matrix4to1[a, b, zeromxn, zeromxn];

        {X, H} = RowHermiteForm[F, oplist];
        {{X11, X12}, {X21, X22}} = matrix1to4[X, {n, n}];

        If [ allflag,
             Print[ "Hermite Form (right side):" ];
             Print[ MatrixForm[H] ];
             If [ ! SameQ[F . X, H],
              Print[ "Error:  F . X != H" ] ] ];

        a . X12 ] /;
    Dimensions[a] == Dimensions[b] && ResamplingMatrix[a] &&
      ResamplingMatrix[b]

(*  NormalizeSamplingMatrix  *)
NormalizeSamplingMatrix[ sampmat_ ] :=
    Block [    {diag},
        diag = Map[ vectorgcd, sampmat ];
        { DiagonalMatrix[diag], mat / diag } ]

(*  RegUnimodularResMatrixQ  *)
RegUnimodularResMatrixQ[mat_] := ResamplingMatrix[mat] && Abs[Det[mat]] == 1

(*  ReorderResampling  *)
ReorderResampling[ m_, nvars_, nlist_ ] :=
    ReorderResampling[ m, nvars, nlist, nvars ]
ReorderResampling[ m_, nvars_, nvars_, zlist_ ] :=
    {m, nvars, zlist}
ReorderResampling[ m_?Matrix, nvars_List, nlist_List, zlist_List ] :=
    Block [    {i, j, locs, mnew, nnew, t, znew},
        locs = Flatten[ Map[Position[nlist, #1]&, nvars] ];
        znew = Part[zlist, locs];
        dim = Length[nvars];
        mnew = m;
        nnew = nvars;
        For [ i = 1, i <= dim, i++,
          For [ j = i+1, j <= dim, j++,
            If [ locs[[i]] > locs[[j]],
             t = locs[[i]]; locs[[i]] = locs[[j]]; locs[[j]] = t;
             t = znew[[i]]; znew[[i]] = znew[[j]]; znew[[j]] = t;
             t = nnew[[i]]; nnew[[i]] = nnew[[j]]; nnew[[j]] = t;
             mnew = permuteMatrix[mnew, i, j] ] ] ];
        {mnew, nnew, znew} ] /;
    SameQ[ Dimensions[m],
           { Length[nlist], Length[nlist] },
           { Length[zlist], Length[zlist] } ]

permuteMatrix[m_, i_, j_] :=
    Block [    {mnew = m, t},
        (* switch rows *)
        t = mnew[[i]]; mnew[[i]] = mnew[[j]]; mnew[[j]] = t;

        (* switch i,j entries in each row *)
        t = mnew[[i,j]]; mnew[[i,j]] = mnew[[i,i]]; mnew[[i,i]] = t;
        t = mnew[[j,j]]; mnew[[j,j]] = mnew[[j,i]]; mnew[[j,i]] = t;

        mnew ]

(*  RelativelyPrime  *)
RelativelyPrime[a_Integer, b_Integer, rest___] :=
    ( GCD[a, b] == 1 )
RelativelyPrime[a_?PolynomialQ, b_?PolynomialQ, rest___] :=
    ( PolynomialGCD[a, b] == 1 )
RelativelyPrime[a_?DiagonalMatrixQ, b_?DiagonalMatrixQ, rest___] :=
    Block [    {i, n, prime = True},
        n = First[ Dimensions[a] ];
        For [ i = 1, i <= n, i++,
              If [ GCD[ a[[i,i]], b[[i,i]] ] != 1,
               prime = False; Break[] ] ];
        prime ] /;
    ResamplingMatrix[a] && ResamplingMatrix[b] &&
      Dimensions[a] == Dimensions[b]
RelativelyPrime[a_?ResamplingMatrix, b_?ResamplingMatrix, Left] :=
    RegUnimodularResMatrixQ[ GCRD[a, b] ] /;  (* or does gcrd = Imat? *)
    Dimensions[a] == Dimensions[b]
RelativelyPrime[a_?ResamplingMatrix, b_?ResamplingMatrix, Right] :=
    RegUnimodularResMatrixQ[ GCLD[a, b] ] /;  (* or does gcld = Imat? *)
    Dimensions[a] == Dimensions[b]

(*  ResamplingMatrix  *)
intcond[a_Integer] := True
intcond[a_?NumberQ] := False

intcond/: Format[ intcond[a_] ] := StringForm[ "Head[``] == Integer", a ]

ResamplingMatrix[a_] :=
    MatrixQ[ a ] && Apply[ SameQ, Dimensions[a] ] &&
    Apply[ And, Map[intcond, Flatten[a]] ] && ( Det[a] != 0 )

(*  ResamplingMatrixMod  *)
ResamplingMatrixMod[n_?IntegerVectorQ, sampmat_?MatrixQ] :=
    ResamplingMatrixMod[n, sampmat, Inverse[sampmat] ] /;
    ( ! SameQ[ ResamplingMatrix[sampmat], False ] ) &&
    Apply[ SameQ, Prepend[Dimensions[sampmat], Length[n]] ]

ResamplingMatrixMod[n_, sampmat_, invsampmat_] :=
    n - ( sampmat . Floor[ invsampmat . n ] )

(*  RowHermiteForm  *)
Options[RowHermiteForm] := { Dialogue -> False }
RowHermiteForm[ mat_, op___ ] :=
    Block [    {Ht, oplist, Xt},
        oplist = toList[op] ~Join~ Options[RowHermiteForm];
        {Xt, Ht} = ColumnHermiteForm[ Transpose[mat], oplist ];
        { Transpose[Xt], Transpose[Ht] } ] /;
    MatrixQ[mat] && IntegerVectorQ[ Flatten[mat] ]

(*  UpsampledSystem  *)
linearmap = {}

addvector[ a_ ] := linearmap = a /; EmptyQ[linearmap]
addvector[ a_ ] := a /; MemberQ[linearmap, a]
addvector[ a_ ] := Block [ {}, AppendTo[linearmap, a]; a ]

getcoeffs[ expr_, vlist_ ] :=
    Block [ {coeffs},
        coeffs = Map[ Coefficient[expr, #]&, vlist ];
        {expr - coeffs . vlist, coeffs} ]

testcoeffs[ {leftover_, coeffs_}, vlist_, zerovector_ ] :=
    ! SameQ[coeffs, zerovector] && MyFreeQ[ {leftover, coeffs}, vlist ]

mDlinearMapping[ expr_, upvarlist_, coeffun_ ] :=
    Block [ {dims, uppair, zerovector},
        dims = Length[upvarlist];
        zerovector = Table[0, {dims}];
        linearmap = {};
        Collect[expr, upvarlist] /.
        ( a_ :> addvector[ coeffun[a, upvarlist][[2]] ] /;
          testcoeffs[coeffun[a, upvarlist], upvarlist, zerovector] );
        linearmap ]

notcoupled[upmatrix_, rowcol_, dims_] :=
    Block [    {oneElementVector},
        oneElementVector = Table[0, {dims}];
        oneElementVector[[rowcol]] = _;
        MatchQ[upmatrix[[rowcol]], oneElementVector] &&
          MatchQ[Transpose[upmatrix][[rowcol]], oneElementVector] ]

UpsampledSystem[ x_, upvarlist_, coeffun_:getcoeffs ] :=
    Block [    {dims, coupledvars, i, j, maxi, pos, rowcol, temp, upmatrix},
        upmatrix = mDlinearMapping[x, upvarlist, coeffun];
        maxi = dims = Length[upvarlist];
        coupledvars = upvarlist;
        rowcol = 1;
        For [ i = 1, i <= maxi, i++,
              If [ notcoupled[upmatrix, rowcol, dims],
               --dims;
               pos = {rowcol, rowcol};
                   upmatrix = Drop[upmatrix, pos];
               upmatrix = Transpose[Drop[Transpose[upmatrix], pos]];
               coupledvars = Drop[coupledvars, pos],
               rowcol++ ] ];
        For [ i = 1, i <= dims, i++,
              If [ SameQ[upmatrix[[i,i]], 0],
               For [ j = 1, j <= dims, j++,
                 If [ ! SameQ[i,j] &&
                      ! SameQ[upmatrix[[j,i]], 0],
                      temp = upmatrix[[i]];
                      upmatrix[[i]] = upmatrix[[j]];
                      upmatrix[[j]] = temp ] ] ] ];
        {upmatrix, coupledvars} ]


(*  S M I T H     D E C O M P O S I T I O N     R O U T I N E S  *)

(* SmithNormalForm *)
Options[ SmithNormalForm ] := { Dialogue -> False }

SmithNormalForm[ mat_?MatrixQ, op___ ] :=
    Block [    {curmat, dialogue, dim, endflag, identm, identn,
         m, n, oplist, r, submatleft, submatright, transmatleft,
         transmatright, uinv, vinv },

        oplist = toList[op] ~Join~ Options[SmithNormalForm];
        dialogue = Replace[Dialogue, oplist];

        curmat = mat;
        {m, n} = Dimensions[mat];
        uinv = identm = IdentityMatrix[m];
        vinv = identn = IdentityMatrix[n];
        r = Min[m, n];

        For [ dim = 0, dim < r - 1, dim++,

              If [ dialogue || SameQ[dialogue, All],
               Print[]; Print[ "Iteration #", dim + 1 ]; Print[] ];

              endflag = False;
              While [ ! endflag,

                  (* part a: move min element to (1,1) *)

                  { minvalue, transmatleft, transmatright } =
                SmithFormA[curmat, dim, m, n, identm, identn];
                  curmat = transmatleft . curmat . transmatright;

                  SmithFormTrace[dialogue, "Part A", curmat];

                  (* part b: eliminate elements on first row/col *)

                  { submatleft, submatright } =
                SmithFormB[ curmat, dim, minvalue,
                        m, n, identm, identn ];
                  curmat = submatleft . curmat . submatright;

                  SmithFormTrace[dialogue, "Part B"];
                  SmithFormTrace2[dialogue, "left", submatleft];
                  SmithFormTrace2[dialogue, "right", submatright];
                  SmithFormTrace[dialogue, " ", curmat];

                  (* update U^-1 and V^-1 *)
                  uinv = submatleft . transmatleft . uinv;
                  vinv = vinv . transmatright . submatright;

                  (* check for the ending condition *)
                  endflag = SmithFormEndTest[curmat, dim, m, n] ] ];

    {uinv, curmat, vinv} ]

(* SmithOrderedForm --  sort the diagonal matrix *)
Options[SmithOrderedForm] := Options[SmithNormalForm]

SmithOrderedForm[ mat_List, op___ ] :=
    Block [    {d, dim, i, identm, identn, j, jend, m, n, r,
         temp, utransmat, vtransmat, uinv, vinv},

        {m, n} = Dimensions[mat];

        (* First, find its Smith Normal Form decomposition *)
        { uinv, d, vinv } = SmithNormalForm[ mat, op ];

        (* Factor out the signs of the diagonal elements *)
        (*   the diagonal values will all be positive    *)
        (*   propagate the sign info to U^-1         *)
        Clear[i];
        uinv = Sign[ Table[ Take[d[[i]], m], {i, 1, m} ] ] . uinv;
        d = Abs[d];

        (* Sort the elements along the diagonal *)
        utransmat = identm = IdentityMatrix[m];
        vtransmat = identn = IdentityMatrix[n];
        r = Min[m, n];
        For [ i = 1, i <= r, i++,
              jend = r - i;
              For [ j = 1, j <= jend, j++,
                If [ d[[j,j]] > d[[j+1,j+1]],
                 temp = d[[j,j]];
                       d[[j,j]] = d[[j+1,j+1]];
                       d[[j+1,j+1]] = temp;

                       temp = utransmat[[j]];
                       utransmat[[j]] = utransmat[[j+1]];
                       utransmat[[j+1]] = temp;

                       temp = vtransmat[[j]];
                       vtransmat[[j]] = vtransmat[[j+1]];
                       vtransmat[[j+1]] = temp ] ] ];

        (* We switched rows instead of columns, so transpose it *)
        vtransmat = Transpose[vtransmat];

        (* We have changed  U D V  to  U Tu^-1 Tu D Tv Tv^-1 V    *)
        (*   Dnew = Tu D Tv  which was done in the sort above    *)
        (*   Unew = U Tu^-1  ==>  Unew^-1 = Tu U^-1        *)
        (*   Vnew = Tv^-1 V  ==>  Vnew^-1 = V^-1 Tv        *)
        uinv = utransmat . uinv;
        vinv = vinv . vtransmat;

        (* return U^-1, D, and V^-1 as a list *)
        { uinv, d, vinv } ]

(* SmithReducedForm *)
(* -- A is the original matrix and A = U D V  *)
(* -- D has only positive values              *)
(* -- start with        U^-1 A V^-1   =   D   *)
(* -- iterate with    G U^-1 A V^-1 H = G D H *)
Options[SmithReducedForm] := Options[SmithNormalForm]

SmithReducedForm[ mat_List, op___ ] :=
    Block [    {adj, d, di, gcd, g, h, i, identm, identn,
         lambda, lastd, lasti, lcm, m, mu, n, r, uinv, vinv},

        (*  Find the dimensions of the matrix *)
        { m, n } = Dimensions[mat];
        r = Min[m, n];
        identm = IdentityMatrix[m];
        identn = IdentityMatrix[n];

        (*  Put the matrix in Smith Ordered Form  *)
        { uinv, d, vinv } = SmithOrderedForm[mat, op];

        (*  Convert the diagonal  *)
        lastd = d[[1,1]];
        lasti = 1;
        For [ i = 2, i <= r, i++,
              di = d[[i,i]];
              (* gcd = GCD[lastd, di];
                 {lambda, mu} = bezoutwithgcd[lastd, di, gcd]; *)
              {gcd, {lambda, mu}} = ExtendedGCD[lastd, di];
              If [ gcd == lastd, lastd = di; Continue[] ];

              g = identm;
              h = identn;
              lcm = LCM[lastd, di];

              (* define g matrix *)
              g[[lasti, lasti]] = 1;
              g[[lasti, i]] = 1;
              g[[i, lasti]] = -mu di / gcd;
              g[[i, i]] = lambda lastd / gcd;

              (* define h matrix *)
              h[[lasti, lasti]] = lambda;
              h[[lasti, i]] = -di / gcd;
              h[[i, lasti]] = mu;
              h[[i, i]] = lastd / gcd;

              (* adjust diagonal elements *)
              d[[lasti, lasti]] = gcd;
              d[[i,i]] = lastd = lcm;

              (* update uinv and vinv *)
              uinv = g . uinv;
              vinv = vinv . h;

              lasti = i ];

        { uinv, d, vinv } ]

(* SmithFormSameV *)
SmithFormSameV[ smithfun_, m1_, m2_, op___] :=
    Block [    {d1, d2, diag, u1inv, u2, u2d2, u2inv, v1inv, v2inv},

        {u1inv, d1, v1inv} = smithfun[ m1, op ];

        v2inv = v1inv;
        u2d2 = m2 . v2inv;            (* U2 D2 = M2 V2^-1 *)
        diag = Map[vectorgcd, Transpose[u2d2]];    (* D2 is the gcd of *)
        d2 = DiagonalMatrix[diag];        (*   cols of U2 D2  *)
        u2 = u2d2 . DiagonalMatrix[1/diag];    (* U2 = U2 D2 D2^-1 *)
        u2inv = Inverse[u2];

        {Abs[Det[u2inv]] == 1, {u1inv, d1, v1inv}, {u2inv, d2, v2inv}} ]

(* SmithFormSameU *)
SmithFormSameU[ smithfun_, m1_, m2_, op___ ] :=
    Block [    {d1, d2, diag, u1inv, u2inv, u2invd2, v1inv, v2inv},

        {u1inv, d1, v1inv} = smithfun[ m1, op ];

        u2inv = u1inv;
        d2v2 = u2inv . m2;            (* D2 V2 = U2^-1 M2 *)
        diag = Map[vectorgcd, d2v2];        (* D2 is the gcd of *)
        d2 = DiagonalMatrix[diag];        (*   rows of D2 V2  *)
        v2 = DiagonalMatrix[1/diag] . d2v2;    (* V2 = D2^-1 D2 V2 *)
        v2inv = Inverse[v2];

        {Abs[Det[v2inv]] == 1, {u1inv, d1, v1inv}, {u2inv, d2, v2inv}} ]


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ]


(*  H E L P     I N F O R M A T I O N  *)

Block [    {newfuns},
    newfuns =
      { BezoutNumbers,        CommutableResamplingMatricesQ,
        DiagonalMatrixQ,        DistinctCosetVectors,
        EuclidFactors,
        GCLD,            GCRD,
        ColumnHermiteForm,        ImpulseTrain,
        IncList,            IntegerVectorQ,
        LCLM,            LCRM,
        NormalizeSamplingMatrix,    RegUnimodularResMatrixQ,
        RelativelyPrime,        ReorderResampling,
        ResamplingMatrix,        ResamplingMatrixMod,
        RowHermiteForm,
        SmithFormSameU,        SmithFormSameV,
        SmithNormalForm,        SmithOrderedForm,
        SmithReducedForm,        UpsampledSystem };
    Combine[ SPfunctions, newfuns ];
    Apply[ Protect, newfuns ] ]


(*  E N D     M E S S A G E  *)

Print["Functions supporting multirate signal manipulations have been loaded."]
Null